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^ ; ABSTRACT 

, Direct detection of low-frequency gravitational waves (10^^ — 10^^ Hz) is the main 

goal of pulsar timing array (PTA) projects. One of the main targets for the PTAs is 
to measure the stochastic background of gravitational waves (GWB) whose charac- 
teristic strain is expected to approximately follow a power-law of the form hc{f) ~ 
A(//yr~^)", where / is the gravitational-wave frequency. In this paper we use the 
current data from the European PTA to determine an upper limit on the GWB am- 
5h ' plitude ^ as a function of the unknown spectral slope a with a Bayesian algorithm, 

by modelling the GWB as a random Gaussian process. For the case a = —2/3, which 
is expected if the GWB is produced by supermassive black-hole binaries, we obtain a 
95% confidence upper limit on ^ of 6 x 10~^^, which is 1.8 times lower than the 95% 
confidence GWB limit obtained by the Parkes PTA in 2006. Our approach to the data 
analysis incorporates the multi-telescope nature of the European PTA and thus can 
serve as a useful template for future intercontinental PTA collaborations. 

Key words: gravitational waves - pulsars: general - methods: data analysis 
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1 INTRODUCTION 

The first direct detection of gravitational waves (GWs) 
would be of great importance to astrophysics and fun- 
damental physics: it would confirm some key predictions 
of general relativity, and lay the foundation for observa- 
tional gravitational-wave astronomy. Pulsar Timing Arrays 
(PTAs) are collaborations which aim to detect low-frequency 

* Email: haasteren@strw.leidenuniv.nl 



(10~^ — 10~*Hz) extragalactic gravitational waves directly, 
by using a set of Gala ctic millisecond pulsar s as nearly- 
perfect Einstein clocks (|Foster fc Backeil Il990l ). The basic 
idea is to exploit the fact that millisecond pulsars create 
pulse trains of exceptional regularity. GWs perturb space- 
time between the pulsars and the Earth, and this creates 
detectable deviations from the stri ct periodicity in the ar- 
rival times of the pulses (TOAs) jEstabrook fc Wahlguistl 
Il975l : ISazhinll 19781 : lDetweileilll979l ). 

One of the main astrophysical targets of the PTAs 
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is to measure the stochastic background of gravitational 
waves (GWB). This GWB is expected to be generated 
by a large number of black-hole binaries located at the 
centres of galaxies (iBegelman et alj 1980l: Phinnev| 200 ll: 
Jaffe fc Backeil l2003l : IWvithe fc Loebl l200a: ISesana et al l 
2008h . by relic gravitational waves (jGrishchukl l2005l l. 



or, more speculatively, by oscillating cosmi c-string loops 
(|Damour fc Vilenkinll200a lOlmez et al.ll2010l ). 

Currently, there are three independent PTA groups: 

(i) the Australian-based programme PPTA, the Parkes Pul- 
sar Ti ming Array, which uses data from th e Parkes tele- 
scope (|Hobbs et al.ll2009l : IVerbiest et al]|2010l ) , and archival 
Arecibo data. 

(ii) the North-American based programme NANOGrav, 
North-American Nanohertz Observatory for Gravitational 
waves, which uses both the Gre en Bank Tel escope (GBT), 
and the Arecibo radio telescope (|jenetll2009l ). 

(iii) and the European programme EPTA, European Pul- 
sar Timing Array, which uses five different radio telescopes: 
the Lovell telescope near Manchester, United Kingdom, the 
Westerbork Synthesis Radio Telescope (WSRT) in the north 
of the Netherlands, the Effelsberg Telescope (EPF) near 
Bonn in Germany, the Nangay Radio Telescope (NRT) near 
Nangay in France, and the Sardinia Radio Telescope (SRT) 
in Sardinia, Italjij. 

It is likely that the first detection of GWs by a PTA will oc- 
cur as a result of a joint effort of all current P TA projects: 
an In ternational Pulsar Timing Array (IPTA; iHobbs et al.l 
l2010t l. This will involve the combination of data from sev- 
eral different telescopes, each of them with its own specific 
hardware elements and software analysis tools. Combining 
data of different observatories is a challenging task, which 
requires extra care when dealing with t he high quality data 
of modern observatories (jjansse nl l2009D . 

In this EPTA paper, we present a methodology on how 
to combine the data from several radio telescopes and use it 
in an optimal way to obtain the information on extragalactic 
gravitational waves. We use the data from three different 
radio telescopes located on the European continent, to place 
a new upper limit on the amplitude of the GWB. As part 
of our analysis, we obtain detailed information about the 
statistical properties of the individual pulse time series. 

The calculation of upper limits on the GWB, based 
on pulsar timing, go as far back as t h e early 1990's 
llStinebring et aT I Il990l : iKaspi et al. I Il994l : iMcHugh et al.l 
1 19961 : lLommenll2002l ). These analyses have been based on 
high quality datasets for single millisecond pulsars. The 
most string e nt up per limits have been obtained recently by 
Ijenet et al. (|2006l ). who have used PPTA data and archival 
Arecibo data for several mil lisecond pul s ars. O ur dataset 
is different from that used bv lJenet" et al.l (|2006l ) since it in- 
cludes only the pulse times of arrival measured by the EPTA 
telescopes, even though some of the pulsars are being timed 
by multiple PTA groups. The Bayesian algorithm we use to 
obtain an upper limit on the GWB is also different from 
the algorithms used by all of the previous studies. Its poten- 
tial advantages include the use of cross correlations between 



^ The SRT is e xpected to become operational in 2011 
llTofani et al.ll2008l) 



TOAs of different pulsars, and the simultaneous constraint 
on both the amplitude and spectral index of the GWB. 

The outline of the paper is as follows. In Section [2] we 
give a brief general overview of pulsar timing observations. 
In Section|3]we detail the observations from all of the EPTA 
telescopes which were used for this paper's analysis. We out- 
line the data analysis procedure in Section 31 after which, in 
Section [5] we present the upper limits on the amplitude of 
the GWB, and also the spectral analysis of the individual 
pulsar noises. Finally, in Section |B] we discuss the astrophys- 
ical implications of our results. 



2 EPTA DATA ANALYSIS 

In this section we present a brief overview of the observa- 
tions, instrumentation and data analysis used at the differ- 
ent EPTA observatories for transforming a series of mea- 
sured pulses to a TO A. 

The complete data reduction process that converts the 
incoming data stream from a radio telescope into one single 
TOA per observation, called "the pipeline" , is optimised by 
hand with much care and is observatory specific. The pro- 
cess can be described in five general steps, shown in Figure 

m 

1) The incoming radio waves are received by the telescope. 

2) The signal is converted from analog to digital, at a 
Nyquist sampled rate. 

3) Data is (coherently) de-dispersed and, if possible. Stokes 
parameters are formed. 

4) The de-dispersed timeseries are folded at the pulsar pe- 
riod, resulting in averaged pulse profiles. Typically a times- 
pan containing several 10^ pulses is used for each TOA. 

5) A cross-correlation with a template pulse profile yields a 
TOA and associated uncertainty (jTavlor 199^ ). 

Individual pulse amplitudes and pulse shapes are highly 
irregul ar, and pulse phases vary significantly from pulse to 
pulse (|Cordes fc Shannonll201ol '). Therefore careful averag- 
ing (folding) has to be performed to obtain a single TOA. 
Furthermore, the interstellar medium (ISM) results in signif- 
icant delays of the arrival time of the pulses over the receiver 
bandwidth. As a large bandwidth is required to reliably de- 
tect a pulse, accounting for the ISM is key for precision 
timing. 

Differences in templates used, e.g. the use of inte- 
grated profiles versus analytic templates, all based on single- 
observatory data, and the difference in definition of the ref- 
erence point in a template will result in offsets between data 
sets generated by different observatories. All extra offsets in 
our data will lead to information loss of other signals like the 
GWB. Therefore, using a common template for each pulsar 
at all observatories is desirable, and will be implemented in 
the near future. 

The realisation of the five steps and therefore their out- 
put (the resulting TOA) might differ among observatories. 
Understanding and accounting for those differences is essen- 
tial for the correct analysis and optimal combining of the 
EPTA data. A more detailed study on this subject is in 
preparation (Janssen et al. 2011). 

The cross-correlation between the folded profile and th e 
template yields an uncertainty of the TOA (|Tavloil[l992l ). 
One would like this uncertainty to be solely due to the ra- 
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Output 

Figure 1. The processing pipeline for pulsar timing, step by step 



diometer noise, i.e. the noise intrinsic to the measurement, 
but in practice the errors sometimes appear to have been 
systematicaUy over- or underestimated. It is a common prac- 
tice, which we foUow here, to allow for an extra parameter 
to multiply these unce rtainties for each pulsar- observatory- 
backend combination (jHobbs fc Edwardal2006l ). This extra 
multiplicative factor allows the TOA uncertainties to sta- 
tistically account for the TOA scatter: the deviations of the 
strict periodicity of the pulses. This is clearly unsatisfactory, 
and in future timing experiments the origin of the predicted 
and measured TOA scatter will have to be thoroughly in- 
vestigated. 



3 EPTA OBSERVATIONS 

3.1 Overview of the observatories 

We have used pulsar timing observations of five radio pul- 
sars, observed with three of the EPTA telescopes, to set a 
limit on the GWB. See TablefT] Fig.[2] and the Appendix 
for an overview of the data sets used and the properties of 
each telescope. Each pulsar was observed on average once 
every month for 30 minutes at each telescope. Although ad- 
ditional observing frequencies are commonly used at WSRT 
and EPF, their respective 1380 and 1400 MHz observing 
bands have the best sensitivity and result in the highest 
precision TOAs. Therefore we have only used observations 
taken at those frequencies at WSRT and EPP for the analy- 
sis presented in this paper. The data were either coherently 
de-dispersed (NRT and EPP) or incoherently de-dispersed 
(WSRT) . The observations were folded and cross-correlated 
with an analytic template (EPP), or a high S/N, observatory 
specific, template (WSRT & NRT), to calc ulate one time- 
of-arr ival (TOA) per observation. See e.g. iLazaridis et al.l 
for a more complete description of the observing pro- 
cedures and data analysis at the different observatories. 

As discussed, any change to the pipeline or to the in- 



put of the pipeline can result in a difference in the calculated 
TOAs. We emphasise that it is essential to correctly identify 
these systematic effects and include them in the modelling 
of the TOAs. In our analysis, we have done this by intro- 
ducing jumps between TOAs of the same pulsar anywhere 
the pipeline differs in some way. 

Once the complete set of data for each pulsar is ob- 
tained, and corrected for global drifts by comparing to UTC, 
it is fit with the timing model. The timing model is a multi- 
parameter fit that represents our best knowledge of the 
many deterministic processes that influence the values of 
the TOAs. The timing residuals are then produced by sub- 
tracting the timing model, which is subsequently optimised 
by minimising these residuals through a least-squares fit. 
This was d one using the pulsar timing package Tempo2 
ijHobbs et al..i2006 ). 

3.2 Selection of data sets 

The European observatories have been timing millisecond 
pulsars for many years, and potentially all of that data could 
be used in th e calculation of an upper limit on the GWB. 
However, like lJenet et al.l (|2006l l we choose to use only the 
data from the pulsars which perform best as ideal clocks, 
e.g. those with the highest precision TOAs and the most 
straight-forward noise characteristics. 

TOA precision is not the only factor that determines the 
sensitivity to the GWB; other factors like the total timing 
baseline and the number of observations (i.e. TOAs) affect 
this sensitivity as well. A great advantage of the EPTA data 
is that several pulsars have been monitored for a relatively 
long time: over 10 years. To determine which timing resid- 
uals (i.e. pulsar-observatory combinations) are most useful 
for GWB detection, we analyse each dataset separately. By 
doing this we can determine the sensitivity to the GWB of 
a set of TOAs: the lower the 3-(t upper limit /i™'"'(lyr) we 
get using only a particular set of TOAs, the more sensitive 
that set of TOAs is to the GWB. 
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Telescope 


WSRT 


NRT 


EPF 


Equivalent dish size (m) 


93.5 


94.4 


100 


Centre observing frequencies (MHz) 


1380 


1398, 2048 


1400 


Observing bandwidth (MHz) 


80 


64/128 


28-112 


Obs. time per month per pulsar 


1x30 min 


4-6x60min 


1x30 min 


Pulsar backend 


PuMal 


BON 


EBPP 


Dedispersion 


incoherent 


coherent 


coherent 


Used templates 


integrated profiles 


integrated profiles 


analytic 



Table 1. Details of the different EPTA observatories relevant for this work. The NRT observing bandwidth has doubled to 128 MHz in 
July 2009. 



Date (Year) 

2000 2002 2004 2006 2008 2010 2012 

I i f J1909-3744 (nrt) 

■I-l ^^■■ii---rJ"»-'fT*P'^- J1744-1 134 (nrt) 

I I TI 
^i#I'-J^^■-•."■■Ii^■■^Vf=_J:^I - 'k- -I- J1744-1134 (etf) 

■I-l j:^'^{][|^i=#*|^-4|5iiifp- i|=---}e^I%t*Ia-isf J1713+0747 (wsrt) 
■I-! jiiis.iri.-..J^J.i. ... ..J_fi. J1713+0747 (etf) 

.1 .! I \- 1 .1 1 ^^B^^&I J1 01 2+5307 (nrt) 



J061 3-0200 (nrt) 




52000 53000 54000 

MJD (Days) 



4 DATA ANALYSIS 

The analysis presente d in this paper broadly follow s the pro- 
cedure introduced in Ivan Haasteren et al.1 (|2009l . vHLML). 
The vHLML Bayesian algorithm relies on creating the 
parametrised models of the timing residuals, and forming 
a probability distribution function (PDF) as a function of 
the model parameters. All known systematic contributions 
of known functional form should be included in the model. 
In the examples used by vHLML the model for the system- 
atic errors included only the quadratic contribution to the 
TOAs from pulsar spindowns. The multi-telescope nature of 
the EPTA requires more complete models for timing resid- 
uals than the one used in vHLML. In this section we show 
how to build and implement these models in practice. 

We first briefly review the algorithm of vHLML in Sec- 
tion 14.11 and 14.21 We then present the extended model we 
use for the analysis of the TOAs in Section [4.31 after which 
we show how we handle TOAs coming from different obser- 
vatories in Section [4.41 



Figure 2. The timing residuals of all the pulsars used in the GWB 
limit calculation. The time in MJD is shown on the x-axis. On 
the left of the dash-dotted line we have placed a sample residual 
with an uncertainty of 1 /^s. 



4.1 Brief review of the vHLML algorithm 

The set of TOAs from all pulsars forms the basic input 
used in the Bayesian data analysis. Many processes influence 
the measured TOA values; in this work we discriminate be- 
tween deterministic processes, like quadratic spindown, and 
stochastic processes, like timing noise: 



The timing residuals of the selected pulsars are shown 
in Figure [2] These flve pulsars significantly outperform the 
other pulsars being timed by the EPTA in terms of how 
well they can limit the GWB amplitude: these five pulsars 
can each individuaUy limit the GWB weU below fec(lyr) — 
10"^'^ for a = —2/3, whereas other current EPTA datasets 
typically perform worse by a factor of several. Since there 
is such a difference between this set of five pulsars, and the 
other pulsars that have been observed by the EPTA, we do 
not expect to gain any significant sensitivity by including 
more pulsars that cannot meet this constraint. We therefore 
choose /i™''''(lyr) < 10"" with a = -2/3 as a constraint 
for including a dataset in our calculation. 

In addition to this constraint, we also demand that 
datasets that just barely satisfy /i™^''(lyr) ^ 10"^"^ do not 
show prominent low- frequency ( "red") timing noise. Our cri- 
terion for presence of the latter is a peak in the posterior 
distribution which is inconsistent with zero amplitude for 
a < 0. 



''(ai) 



,dGt I r.stoch 



(1) 



where represents the i-th TOA of pulsar a, is the 
corresponding contribution to the TOA solely due to de- 
terministic processes, and is the contribution due to 
stochastic processes. 

The effects of deterministic processes are described by 
the set of model parameters ff: t'^^l-^ = is done in 

vHLML, we assume that the stochastic processes are Gaus- 
sian, though their spectra are not necessarily white. In such 
a model, the stochastic processes can be represented by the 
correlation matrix 



/ c I stocli c I stoclix 
\<J^(ai) "^(6i) / 



= C7, 



(ai)(bj) 



= C7, 



(a, 



(2) 



where ^ are the model parameters. 

The key distribution used in a Bayesian analysis is 
the likelihood function, the probability distribution of the 
data for a given model and its parameters. As described in 
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vHLML, for PTAs the likelihood takes the following form: 



4.3 Used model for the TOAs 



P{5t 



^(27r)"detC 



(3) 



exp 



(ai)(i>j) 



where 6^ = (77, ^) , and (5t is the difference between the ob- 
served TOAs, and the fitted TOAs. A Bayesian analysis 
assigns prior distributions Po{0) to the model parameters, 
and explores the parameter space of the posterior distri- 
bution (short-handed simply as the posterior): P{9 \ 5t) — 

L{e)Po{e). 



4.2 Obtaining a marginalised posterior 
distribution 

The posterior P{6 \ St) contains information about all model 
parameters. We need to express the posterior as a function of 
only those parameters that represent the GWB. This process 
is called marginalisation, and consists of integrating over all 
other parameters. The resulting marginalised posterior is the 
posterior probability density of the GWB parameters. 

Marginalisation of a posterior in a high-dimensional 
parameter space is non-trivial, and a direct numerical in- 
tegration is prohibitively computationally expensive. As 
in vHLML, we employ a mix of analytic integration and 
Markov Chain Monte Carlo (MCMC) methods to accom- 
plish this. The marginalisation remains the computational 
bottleneck for the method's effectiveness, as the computa- 
tional time scales with n"^, with n the total number of TOAs 
to be analysed. 

A computational shortcut can be used by analytically 
marginalising over the parameters of the timing model. As 
shown in vHLML, this is possible provided that the parame- 
ters represent signals of known functional form. This condi- 
tion is equivalent to the requirement that the timing residu- 
als generated by the timing model are linear with respect to 
its parameters: St — d{a — a), where St is the timing resid- 
ual, d is a proportionality constant, a is the best fit value for 
the model parameter, and a is the model parameter. While 
this is always true for quadratic spindown as considered ex- 
plicitly in vHLML, it is generally not true for other timing 
model parameters. However, when the deviations of the tim- 
ing model parameters from their best-fit values are small, it 
is a good approximation that the residuals generated by the 
timing model are linear with respect to the deviations from 
their best-fit values: St ~ d{a — a). 

Analytically marginalising over the timing model is 
therefore possible, and by doing so the number of param- 
eters that must be integrated over numerically by the use 
of MCMC is reduced greatly. Dependent on the model we 
use to describe the statistics of the timing residuals, the 
number of parameters left to explore is then just several per 
pulsar/backend combination. The results of the analysis can 
be presented as a marginalised posterior as a function of any 
parameter in the model, provided that this parameter was 
present in the MCMC run. 



We divide the actual parameterisation in 3 parts: 

a) The deterministic timing model. 

b) The gravitational-wave background. 

c) Other stochastic processes (e.g., timing noise). 

In this section we discuss how we have taken these into ac- 
count in our data analysis. 

As a first step, the TOAs are processed using the soft- 
ware package Tempo2, in order to determine the best-fit 
timing model. This procedure consists of the following steps: 

1. Tempo2 requires an initial guess aoi for the timing model 
parameters a; in order to find timing residuals (pre-fit tim- 
ing residuals). 

2. It then constructs an approximation to the timing model, 
in which the timing residuals depend linearly on ai — aoi- 

3. It finds the best-fit at within this linear approximation, 
and uses those values to update the timing residuals using 
the full non linear timing model (post-fit timing residuals). 

4. The newly obtained parameters and corresponding tim- 
ing residuals are then judged by the person performing the 
model fitting, and if determined necessary the newly ob- 
tained parameters can act as the initial guess for a new 
fitting iteration. Tempo2 also allows adjustment and fitting 
of Oi one by one. 

Finding the timing solution with Tempo2 is not fully al- 
gorithmic, but typically requires someone experienced with 
pulsar timing analysis, who approaches the TOAs fitting in 
several different ways, which ensures that phase coherence 
is maintained and that the relevant deterministic model pa- 
rameters are included properly. Though this strategy works 
well in practice, we should remain conscious of the possi- 
bility that different solutions might be obtained by differ- 
ent observers, who may also choose to include additional 
model parameters. In the appendix we present the tim- 
ing solutions we found for the analysed pulsars. These are 
the values we used as our initial guess, aoi- Note that these 
aoi and their uncertainties, although created with Tempo2 
using the same datasets that we base our upper limit on, 
do not include our model for the red noise. The values and 
uncertainties we list in the appendix therefore do not repre- 
sent our best estimates if we were to take into account the 
red timing noise. Although calculating these best estimates 
of Qi is reasonably straightforward, these estimates are not 
accessible in our MCMC because we have marginalised over 
these parameters analytically. The calculated upper limit on 
the GWB, however, does include all these effects, and there- 
fore automatically incorporates the removal of power from 
the low-frequency GW signal by fitting for the timing model 
parameters and jumps. 

In the above mentioned step 2 where the timing model 
is linearised, we have made an important simplification that 
we now describe in more detail. Since we take into account, 
and marginalise over, all timing model parameters in our 
algorithm, we are effectively working with the TOAs instead 



^ Qualitatively, experienced observers are rightfully so very con- 
fident in their timing solutions. Quantitatively however, the only 
statistical tool currently available for observers to check whether 
the timing solution is reasonable is the reduced statistic. But 
since the error bars obtained with the cross-correlation technique 
cannot be fully trusted, the same holds for the statistic. 
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of just the timing residuals. However, the timing model has 
been linearised by Tempo2 with respect to at — aoi- This 
implies that we need to be sufficiently close to aoi in the 
parameter space for this approximation to be valid, which 
means that the timing residuals derived with Tempo2 need 
to be approved by the person fitting the data, before using 
these as inputs in the Bayesian algorithm. 

The stochastic component contributing to the TOAs is 
characterised as follows. Firstly, general relativity describes 
how the timing residuals of a pair of pulsars are correlated 
due to gravitational waves: 

3 1 - cos 6la6 , / 1 - cos 6'a(, \ 1 1 - COS 6*06 1 1 

U - 2 2 2 )~l 2 +2 + 2^"'" 

(4) 

where Ogfj is the angle bet ween pulsar a and pul- 
sar b ijHeUings fc DownsI Il983l l. The G WB spectrum is 



where C°^i)(6j), C{ai){bj) 



are the correlation 



parametrised a s a power-law of the form iMaggiore 



Sesana et aLll2008l ) 



PhinnevI [iooiyjaife fc Backeij l2003l: IWvithe fc Loeb 



200c 



200a 



A 



(5) 



were he is the characteristic strain as used in Ijenet et al.l 
()2006f l. A is the amplitude of the signal, and a is the spectral 
index. This then results in a correlation matrix for the GWB 
(vHLML): 



GW 
(ai)(f)j) 



(2-)^f^ 

00 



{V{~2 + 2a) cos (yra) [fLrf 



Ult) 



(2n)! (2n + 2q - 2) 



(6) 



where, as in vHLML, r = 2Ti\ti — tj\, and /_l is a cut-off 
frequency, set much lower than the lowest GW frequency we 
are sensitive to. 

Secondly, the stochastic timing noise for each individual 
pulsar is split into three components: 

1) Individual errors of TOA determination from the cross- 
correlation, represented by the TOA error bars. An extra 
free parameter, called the EFAC value, is commonly intro- 
duced by pulsar observers in order to a ccount for possible 
mis-c alibration of the radiometer noise (|Hobbs fc Edwarda 
I2OO6I ): this parameter is a multiplier for all of the TOA error 
bars for a given pulsar. 

2) An extra white noise component, independent of the error 
bars. This basically acts as extra non-time-dependent noise, 
and the parameter is often called an EQUAD parameter. 

3) Red noise, consisting of a power-law spectrum in the tim- 
ing residuals. This component allows for structure in the 
timing residuals. 

All three timing noise components are uncorrelated between 
the pulsars. 

The resulting correlation matrices from components 1, 
2, and 3, as derived in vHLML, are given by: 



{ai){bj) 

,WN _ ^r2 c c 

{ai){bj) — i^aOabOij 

— RaSab 



(ai)(bj) 



(27!-) .Il 



{r(-2 + 2q„) cos (Tvaa) ifLrf-^" 



(2n)!(2n + 2Q-2) 



(7) 



matrices corresponding to the error bars, the extra white 
noise, and the red noise respectively, with a and b denoting 
the pulsar number, i and j denote the observation number. 
At is the TOA uncertainty (the error bar) as calculated in 
the pipeline, Ea is the scaling parameter of the error bars for 
the a'th pulsar (the EFAC factor), Na is the amplitude of 
the white noise, Ra is the amplitude of the red timing noise, 
Qa is the spectral index of the red noise spectrum of pulsar 
a, and t is the time difference between two observations. 



4.4 Combining datasets 

The previous section gives a complete description of the 
model we use to analyse the TOAs of a single pulsar, ob- 
served with one telescope. That model does not yet account 
for the use of different observatories. In this section we ex- 
plain what we do to accomplish this. 

As discussed in Section[2] the reduced data products are 
(sometimes subtlety) influenced by many different compo- 
nents of the reduction process. In order to account for slight 
offsets between TOAs, introduced by using slightly different 
reduction procedures on individual datasets, a calibration 
term needs to be introduced when merging TOAs from dif- 
ferent observing systems. This extra calibration term takes 
the form of a "jump", an arbitrary phase offset between 
datasets, which is fit for simultaneously with other timing 
model parameters. We use the term dataset for any series of 
TOAs that can be analysed without a jump. In practice this 
is any series of TOAs, of the same pulsar, observed with the 
same hardware elements, and processed with the same algo- 
rithms, at the same observing frequency. Here we combine 
7 such datasets (those shown in Figure [2]). 

Jumps have been used routinely when combining 
data of differ ent observatories and/or data recorders (e.g., 
IJansse nl l2009l '). This allows us to find a single solution for the 
timing model of a pulsar timed by multiple observatories. 
However, the TOAs produced by pipelines at different ob- 
servatories may have different statistical properties. In order 
to account for this, we allow the stochastic contributions in 
our model discussed in Section r4.3l to vary between datasets: 

1) One timing model per pulsar (taken directly from 

TEMP02) 

2) Jumps between different datasets 

3) A scaling factor for the error bars (EFAC) for each dataset 

4) An extra white noise component (EQUAD) for each 
dataset 

5) Power law red noise for each dataset 

A major advantage of this approach is that it allows one to 
detect statistical differences between observatories that may 
be introduced by different algorithms/components, and then 
use this feedback to iteratively improve our datasets. 

The analysis of the TOAs consists of two steps. In the 
first step Tempo2 is used to find the timing solution for a 
single pulsar. This includes possible jumps between datasets. 
Once the timing solution is obtained, the results are passed 
on to the Bayesian algorithm. The Bayesian algorithm then 
analytically marginalises all parameters of the timing model, 
including jumps, while using MCMC to explore the rest of 
the parameter space. 
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Now that we have developed the necessary framework to 
analyse the TOAs, we apply the algorithm to the observa- 
tions. In the following sub-sections we explain in detail how 
we selected the five pulsars that we already mentioned in 
Section 13.21 and we present the GWB upper limit we are 
able to calculate using observations of those pulsars. 



5.1 Selecting the most constraining datasets 

For any pulsar, obtaining the timing solution and timing 
residuals is the first step after obtaining the TOAs. The 
timing residuals of the pulsars used in this work are shown in 
Figure [21 and the parameters of the timing model are shown 
in the Appendix. The timing model also includes several 
jumps as some of these pulsars have been observed with 
several European telescopes. The timing solutions we find 
are quite consistent with the values already published in the 
literature. Given that we are solving for 56 parameters, it 
is to be expected that one or two parameters deviate at the 
2-0" level. The only unexpected outlier we find is the proper 
moti on in right as c ension of J1713-I-0747, which deviates 
from ISplaver et al.l (|2005l ) by over S-cr. Given that we are 
combining data of several telescopes, and that we do not 
take into account our red noise models in listing these timing 
solutions, we postpone exploring this difference to future 
work where the focus lies on investigating the statistics of 
the timing model parameters in the presence of red noise. 
Such an investigation is beyond the scope of this manuscript. 

With the model of the systematic contributions in place, 
we first perform the analysis separately for each of the 
datasets and obtain the posterior probability distribution 
for their intrinsic noise parameters, specified in Equation 
((TJ of the previous section. Note that at this stage of the 
analysis the contribution from a GWB is not yet included. 
We determine a marginalised posterior for each pulsar as a 
function of the following parameter combinations: 

1) EFAC vs. EQUAD 

2) Red noise amplitude vs. red noise spectral index 

In both cases, the posterior is marginalised over all parame- 
ters but two, and the resulting 2-dimensional distribution is 
displayed as contours at the 1-, 2-, and 3-a level (the regions 
where respectively 68%, 95%, and 99.7% of the volume of 
the posterior is enclosed). 

As an example we consider the TOAs of pulsar 
J1713-I-0747, which consist of data taken with Effelsberg and 
Westerbork. Here we focus on the marginalised posterior 
distributions that represent information about the Effels- 
berg TOAs; these distributions and the residuals are shown 
in Figs [3] and 21 A traditional non-Bayesian analysis of the 
Effelsberg TOAs consists of a fit to the timing model with 
Tempo2, which yields the best-fit parameters, the corre- 
sponding uncertainties, and a reduced statistic. The re- 
duced is defined as: 



2 



(8) 



51200 51900 52600 53300 54000 54700 
lulJD (Days) 




Tempo2 x 
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Figure 3. The marginalised posterior of J1713+0747 (Effels- 
berg), as a function of the EFAC and EQUAD parameters. The 
contours are at the 1, 2, and 3-(t level, indicating a respective 
volume inside that region of 68%, 95%, and 99.7%. 



tainty of 1°"°^ , and e is the EFAC value. It is usual practice to 
set the EFAC such that the reduced = 1, which is accom- 
plished by: e = ^Jx^{e = 1). For the J1713-H0747 Effelsberg 
TOAs, we have x = 1) = 18-9 and therefore e = 4.35. 

As can be seen from Figure (31 a non-zero red noise com- 
ponent is required to describe the TOAs. The EQUAD pa- 
rameter is consistent with 0-amplitude according to Figure 
O while the EFAC parameter is significantly lower than what 



where n is the number of observations, m is the number of 
free parameters in the least-squares fit, tf"^ is the observed 
TOA, tf is best-fit value of the TOA, Oi is the TOA uncer- 



a Tempo2 vx^ estimate would give. This tells us that no 
separate white-noise component is required to describe the 
TOAs: all the uncorrelated scatter can be assigned to the 
error bars on the TOAs. It is also of interest that in this 
case the EFAC parameter is much smaller, and indeed much 
closer to 1, than the more traditional estimator \/x^. The 
data is reasonably well-modelled by just the presence of red 
noise. 

It is also worth noting that, due to practicalities having 
to do with hardware changes at the observatories, the TOAs 
of J1713-I-0747 end at an earlier epoch than the other 4 
pulsars. Although in the future the inclusion of this data 
will obviously benefit the sensitivity to the GWB, we note 
that the GWB limit is not negatively effected by this lack 
of overlap of the TOAs between pulsars. 

The analysis of the TOAs of the other pulsars yields 
similar, but slightly different results. As can be seen in the 
appendix, some pulsars do have non-negligible white noise, 
and some do appear to have EFAC values significantly differ- 
ent from 1. As of yet we do not have a complete explanation 
for the exact form of the marginalised posteriors. 

We present the marginalised posterior as a function of 
the red noise parameters in an intuitive way: as pointed out 
in Section 13.21 we use the same units for the red noise am- 
plitude and red noise spectral index as we use for the GWB 
parameters. For the analysis of TOAs of just one pulsar, the 
red noise can now be thought of as if it was generated solely 
by a GWB with a certain amplitude and spectral index. In 
this case, the marginalised posterior for the red noise pa- 
rameters shows us what upper limit we are able to place on 
the GWB amplitude as a function of spectral index. 

We choose a 3-o- threshold of Ra ^ 10"^'' at a spectral 
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Figure 4. The marginalised posterior of J1713+0747 (Effels- 
berg), as a function of the power-law red noise parameters: the 
amplitude and the spectral index. The contours are at the 1, 2, 
and 3-(T level, indicating a respective volume inside that region of 
68%, 95%, and 99.7%. 



index of = —2/3. Based on the marginalised posteriors 
of all the EPTA pulsars, we can decide whether a particular 
dataset can put a constraint on the GWB lower than this 
or not. Using this threshold we include five pulsars in our 
final analysis. These five significantly outperform the other 
pulsars in terms of how well they can limit the GWB ampli- 
tude, and we do not expect to gain any significant sensitivity 
by including more pulsars in our current archival data sets. 
The residuals of the pulsars we use in our combined analysis 
are shown in Figure (2] More datasets will be added after 
some extensive and detailed recalibration procedure of ex- 
isting datasets. 



5.2 GWB upper limit 

Now that we have selected our pulsars that can significantly 
contribute to a GWB limit, we are in the position to infer 
the amplitude and spectral index of the GWB. Our model of 
the combined data of the five pulsars we selected in Section 
15.11 consists of all sources we included in the analysis for 
the single pulsars, and an extra source that corresponds to 
the GWB. As discussed in Section 14.31 the GWB source 
is a power-law correlated between pulsars as described by 
Equation Q. 

As before, we use MCMC to sample the posterior dis- 
tribution while analytically marginalising over the timing 
model; now the analytic marginalisation happens simulta- 
neously for the timing models of the five pulsars. In Figure 
[5]we present the posterior, marginalised over all parameters 
except the GWB amplitude and spectral index. In the same 
figure we also show the PPTA published values of the GWB 
limit (jjenet et al.ll2006h . For the expected spectral index for 
a GWB generated by a large number of supermassive black- 
hole binaries, a = —2/3, we find a 95% confidence GWB 
upper limit of hc{lyr) ^ 6 x 10"^''. This is smaller by a 
factor of 1.8 than the previously published PPTA limit. 

As a cross-check with other codes, and to verify that 
we are definitely sensitive to the level of the limit we have 
calculated here, we perform an additional test. We use the 



Jenet et al. (2006) 



van Haasteren et ai. (201 1) 



Expected a for GWB from SlylBHBs 



-0.5 
GWB index a 
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Figure 5. The marginalised posterior distribution as a function of 
the GWB amplitude, and spectral index. The contours marked by 
'van Haasteren et al. (2011)' are the results of this work at the l-a 
and 2-(T level, indicating a respective volume inside that region of 
68%, and 95%. The vertical dash-dotted line at a = —2/3 shows 
where we expect a GWB generated by supermassive black-hole 
binaries. The most recent published limits are shown as the three 
upper limit arrows pointing down, marked by 'Jenet et al. (2006)'. 



TEMP02 plug-in GWbkgrd (|Hobbs et all 120091 ') to generate 
simulated timing residuals as produced by a GWB with an 
amplitude of /ic(lyr). We then create a new set of TO As, 
consisting of the values of the simulated timing residuals 
added to the values of the observed TOAs of the five pul- 
sars that we have analysed in this section. We then redo 
the whole analysis. Current PTAs aim to reach sensitivities 
in th e order of hc{lyr) = 10~^® in the future (|jenet et al.l 
l2005l ). which is over five times more sensitive than the limit 
we achieve here. In the case that the GWB just happens 
to be at the 2-cr level of our current limit, we demonstrate 
what such a fivefold increase in sensitivity could do for our 
ability to measure the GWB parameters by adding a signal 
of hc{lyr) = 30 X 10"^^ to our current TOAs. The result is 
shown in Figure [S] We find that the results are consistent 
with the input parameters of the simulated GWB, and that 
we can reliably detect a GWB in this cas(jf]. The values of 
the GWB parameters we have used to simulate the GWB 
lie within the l-cr credible region of Figure [6] 



6 IMPLICATIONS 

The analysis performed in this work puts an upper limit 
on a GWB with a power-law characteristic strain spectrum 
he = j4(//yr~^)". In the literature, upper limits are typi- 
cally quoted for various values of a, where the considered 
Q depends on the physics responsible for generation of the 
GWB. A useful feature of our approach is that we are able 
to measure a for a strong enough GWB (see vHLML for 
a discussion). The extra degree of freedom in our model, 
Q, necessarily changes the interpretation of the posterior to 



We note that, although such a detection is consistent with a 
GWB, we would need more pulsars to exclude the possibility that 
some other effect is causing the correlated signal we detect here. 
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Figure 6. Same marginalised posterior distribution as inO but 
liere we have injected the residuals of a simulated GWB with 
amplitude /ic(lyr) = 30 X 10^^^ in the data. 



the uncertainties of the key model parameters for different 
scenarios, and come up with feiyr ~ 2 x 10"^® — 4 x 10~^^. 
We are less than a factor of two away from this range, so 
we foresee that we can start to rule out some models in the 

near future. 

Two more results of lSesana et all (|2008D are interesting 
with respect to the limit presented in this work. The first is 
that the frequency dependence of the GWB is expected to be 
steeper than a power-law oc f~^^^ for frequencies / > 10~^ 
Hz. The steepness depends on the chosen model. We have 
incorporated a varying spectral index a in our current anal- 
ysis, and since we are not yet able to detect the GWB, we 
postpone a more thorough investigation of the exact depen- 
dence of he on / to later work with even better datasets. 
The second interesting result is that in the frequency range 
of 10"* Hz ^ / ^ 10"^ Hz, the GWB might be dominated 
by single sources. In that case, a search for just a certain 
characteristic strain spectrum is not appropriate, and we 
note that further investigation is required in this regard. 



some extent. We interpret the 2-cr contour in our plot of the 
marginalised posterior as the upper limit on the GWB as a 
function of a. Fixing a and re-evaluating the 2-(t limit based 
on the posterior for A only does not significantly alter our 
results. 

In this section, we briefiy discuss the implications of the 
new upper limits in the context of two different mechanisms 
for generation of the GWB: binaries of supermassive black 
holes, and cosmic strings. 



6.1 Supermassive black hole binaries 

Several authors discuss the characteristic strain spectrum 
generated by an ensemble of supermassive bl ack holes (SMB- 



HBs) distributed throughout the Universe 
1980 : iPhinnevHiooH : I Jaffe fc Backer! l2003l : 



Begelman et al 



Wvithe fc Loeb 



20031 ) ■ They show that the characteristic strain spectrum 
generated by such black hole binaries can well be approxi- 
mated by a power-law: 



he 



h 



lyr 



-2/3 



(9) 



where /iiyr is a model-dependent constant. Though the form 
of the characteristic strain, the power-law, is quite general 
among the different SMBHB assembly models the authors 
use in their work, the parameterisations and assumptions 
about other physical quantities differ between all investiga- 
tors. The predicted /iiyr therefore differs depending on what 
SMBHB asse mbly scenarios were assumed. 

Recently, ISesana et al.l (|2008l ) have extensively investi- 
gated a wide variety of as s embly scenarios, including those 
considered in Ijenet et al.l (|2006l V For our purposes in this 
work, their most important result is an estimate of h\y^ for 
all modelfl In calculating this value, they take into account 



* The model for the GWB that lSesana et al.l ||2008| ) use is a bro- 
ken power-law. Their hiyr therefore has a slightly different mean- 
ing, and our quoted value should be taken as a crude approxima- 
tion. 



6.2 Cosmic strings 

Several authors have suggested that oscillatin g cosmic 
string loops will produce g ravita ti onal waves (|Vilenkinl 



19811: iDamour fc Vilenkiiil l2005l : lOlmez et all |2010| ). 
Damour fc VilenkinI (120051 ) have used a sem -analytical 
approach to derive the characteristic strain he of the GWB 
generated by cosmic strings: 



14 1/2 
C T 



he{f) = 1.6 X 10 

x(V0.65)^'"' 



-1/2 -1/6 



10-6 



off 
1/3 



-7/6 



(10) 



where ^ is the string tension, G is Newton's constant, c is the 
average number of cusps per loop oscillation, p is the recon- 
nection probability, ecfi is the loop length scale factor, and h 
is the Hubble constant in units of 100km s-^Mpc"^. Usually, 
the dimensionless combination G/i is used to characterise the 
string tension. Theoretical predictions of string tensions are 
10"" ^Gfi^ 10"^ (Dam our fc Vilcn kin 2005). 

From the above expression for the characteristic strain 
generated by cosmic strings, we see that this is again a 
power-law, but now with a = —7/6. Using a standard model 
assumption that c = 1, the facts that p and e^a are less than 
one, and that h is expected to be greater than 0.65, we 
can safely use our derived upper limit on he for a — —7/6 
to limit the string tension: G/x ^ 4.0 x 10^^. This already 
places interesting constraints on the theoretical models, and 
in a few years the EPTA will be able to place much tighter 
restrictions in the case of a non-detection of a GWB: with 
only a factor of five decrease of the upper limit, we would 
be able to co mpletely exclude the 1 0"^^ ^ Gfi ^ 10"^ range 
mentioned in iDamouj fc VilenkinI (|2005l ). 



7 CONCLUSION AND DISCUSSION 

In this paper we have developed the methodology on how 
to handle combined PTA datasets of several telescopes and 
how to robustly calculate a corresponding upper limit on the 
GWB. Our Bayesian approach has handled in a straightfor- 
ward way different data sets of varying duration, regularity, 
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and quality. The current upper limit on the GWB, calcu- 
lated with EPTA data, is /ic < 6 x 10"^'' in the case of 
a = —2/3, as predicted for a GWB created by an ensemble 
of supermassive BH binaries. More generally, the analysis 
has resulted in a marginalised posterior as a function of the 
parameters of the GWB: the GWB amplitude and the spec- 
tral index. 

Due to hardware and software upgrades at the EPTA 
observatories, and due to the ever increasing time baseline 
of the data, we expect the sensitivity to increase greatly over 
the next few years. Especially the combination of the EPTA 
data sets with the data of the other PTA projects seems 
promising. 

The raw telescope data must first undergo careful reduc- 
tion and modelling before it can be analysed by the Bayesian 
algorithm. We have provided some discussion of these pro- 
cesses and have motivated our choice of model for the TO As. 
As part of our analysis, we have studied the probability dis- 
tribution of the pulsar noise parameters, and highlighted 
the crucial importance of precise characterisation of the red 
component of pulsar timing noise. 
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APPENDIX A 

Here we show the timing solutions of all datasets used in 
this work, combined with the posterior distributions for the 
timing noise. 



EPTA GWB Limit 11 



0.5 







cd 



0.5 - 








1744 eff 



Tempos efac 



0.5 







1909 nrt 



I — Tempos efae- 



0.5 







_l I I I L. 



1744 nrt 



' — Tempos efac 




la 
2a 
3a 



1 



3 

Efac 



Figure 7. The marginalised posteriors of all datascts, as a function of the EFAC and EQUAD parameters. The contours are at the 1, 
2, and 3-cr level, indicating a respective volume inside that region of 68%, 95%, and 99.7%. For the J1713-0747 posterior, the Tempo2 
estimate is not shown because it has the off-scale value of 4.4. 
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Figure 8. The marginalised posterior of all datasets, as a function of the power-law red noise parameters: the amplitude and the spectral 
index. The contours are at the 1, 2, and 3-it level, indicating a respective volume inside that region of 68%, 95%, and 99.7%. The more 
negative the value of o, the steeper the power-law spectrum, with the spectrum approaching a white spectrum at the right of the plot. 
We also note that the amplitude of the red noise cannot be trivially scaled linearly to an rms value of the timing residuals. 
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NRT 
53443 - 55030 
107 
769 
1.00 
54236 



Pulsar name J0613-0200 J1012+5307 

Fit and data set 

Telescopes used NRT 

MJD range 53367 - 55012 

Number of TOAs 280 

Rms timing residual (ns) 912 

Reduced value 1.00 

Epoch 54189 

Measured Quantities 

Right ascension, a (J2000) 06:13:43.97385(4) 

Declination, S (J2000) -02:00:47.0720(12) 

Pulse freq., u (s"!) 326.600562095168(13) 

Derivative of pulse freq., u {s'^) . . -1.02281(3) X 10"^^ 

PM in RA, fia (masyr-i) 1.90(4) 

PM in DEC, fis (masyr-^) -10.31(9) 

Parallax, tt (mas) — 

Dispersion measure, DM (cm-3pc) 38.77700 

Binary model DD 

Orbital period, Pt (d) 1.19851257534(5) 

Derivative of orbital period, Pt ■ ■ ■ — 

Epoch of periastron, Tq (MJD) . . . 54189.019(6) 

Projected sm. axis of orbit, x (It-s) 1.09144417(8) 

Longitude of periastron, luq (deg) . 47.1(1.6) 

Orbital eccentricity, e 5.47(15) xlO^^ 

Time of ascending node (MJD) ... — 

EPSl (ei), esino; — 

EPS2 (£2), ecoso; — 

Sine of inclination angle, sin i — — 

Companion mass, A'lc {^q) — — 

Assumptions 

Clock correction procedure TT(TAI) 

Solar system ephemeris model .... DE405 



J1713+0747 



EPF & WSRT 
51426 - 54637 

195 

396 

1.13 
53031 



10:12:33.43241(10) 17:13:49.530782(3) 
+53:07:02.665(2) +07:47:37.52343(8) 
190.26783448248(14) 218.811840486637(30) 



-6.1998(4) X 10-1'' 
3.17(7) 
-24.96(9) 

9.0176 

ELLl 
0.60462272322(4) 



0.58181742(13) 



54236.2078302(3) 
1.18(5)xl0-5 
2.20(5) X 10-5 



-4.0836(2) x 10-1'' 
5.017(12) 
-3.96(3) 
0.915(7) 
15.9907 

DD 

67.8253309255(20) 

53014.9592(7) 
32.34242015(7) 
176.2109(12) 
7.49312(13) xlO-5 



Table 2. The timing solutions for the pulsars used in this paper before applying the Bayesian algorithm. These solutions are determined 
using Tempo2, which uses the International Celestial Reference System and Barycentric Coordinate Time. As a result this tim ing model 
must be modified before being used with an observing system that inputs Tempo format parameters. See lHobbs et al" I liooi) for more 
information. Note that the figures in parentheses are the nominal l-cr Tempo2 uncertainties, with EFACs included, and therefore do 
not include the red noise model. In the GWB limit calculation these respective parameters are marginalised over. Also, the dispersion 
measure quoted here results from combining these observations with EPTA data of other frequencies. These DM values are used in 
the dedispersion, but we didn't include all observations in our GWB analysis. We therefore have not fit for the DM here, and an error 
estimate cannot bo given. 
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Pulsar name 


J1744-1134 


J1909-3744 


Fit and data set 


Telescopes used 


EPF & NRT 


NRT 


MJD range 


51239 - 55001 


53366 - 55127 


Number of TOAs 


159 


113 


Rms timing residual (ns) 


444 


134 


Reduced value 


1.05 


1.00 


Epoch 


53120 


54247 


Measured Quantities 


Right ascension, a (J2000) 


17:44:29.391592(7) 


19:09:47.437982(5) 


Declination, <5 (J2000) 


-11:34:54.5762(6) 


-37:44:14.3176(2) 


Pulse freq., u [s'^) 245.426119777227(4) 


339.31568732355(1) 


Derivative of pulse freq. , ;> (s ^ ) . . 


-5.3817(4) X 10-1*5 


-1.614853(8)xl0-i5 


PM in RA, /ia (masyr~^) 


18.817(10) 


-9.490(11) 


PM in DEC, ^ls (masyr-^) 


-9.30(6) 


-35.89(4) 


Parallax, tt (mas) 


2.602(10) 


1.01(7) 


Dispersion measure, DM (cm~'^pc) 


3.13632 


10.37877 


Binary model 




ELLl 


Orbital period, Pj, (d) 




1.53349947490(6) 


Derivative of orbital period. Pi, ■ ■ ■ 




3.5(5)xl0-" 


Epoch of pcriastron, Tq (MJD) . . . 






Projected sm. axis of orbit, x (It-s) 




1.89799108(11) 


Longitude of pcriastron, luq (deg) . 






Orbital eccentricity, e 






Time of ascending node (MJD) . . . 




54247.169903748(15) 


EPSl (ei), esino; 




6.4(5. 5)xl0-* 


EPS2 (£2), ecoso; 




-3(3) X 10-* 


Sine of inclination angle, sin i 




0.9980(3) 


Companion mass, A'lc (^0) 




0.208(7) 


Assumptions 


Clock correction procedure 


TT(TAI) 


Solar system ephemeris model .... 


DE405 



Table 3. Same as table 2. 



